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ABSTRACT 

This report describes a method for calculating minimum 
and constant injection energy, C 3 , values for conic inter- 
planetary trajectories. The planets are assumed to move 
on mutually inclined elliptical orbits about the sun. 
Lambert's theorem is used to solve the time -constraint on 
the two-point boundary value problem, and an iterative 
procedure, with arrival date’ as the isolation parameter, 
is used to determine the desired C 3 . 



Typical results ?are presented for Earth to Jupiter 
trajectories in the 1970-71 launch period. The comparisons 
made with previously used methods indicate the procedures 
described in the report produce results of comparable 
accuracy, with a savings in computation time. J Jit-* 
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CONIC SOLUTIONS TO THE INTERPLANETARY TRANSFER PROBLEM 
WITH CONSTANT OR MINIMUM" INJECTION ENERGY 


SUMMARY 

The procedures for calculating constant and minimum 
injection energy, C 3 , solutions to the time -constrained 
two-point boundary value problem ’of interplanetary flight 
are presented. The boundary values, position and velocity 
of the launch pianet at launch time and target planet at 
arrival time, are determined using mean elliptical elements 
to describe the path of the planets. The time-constraint 
on the problem is solved using a modified form of Lambert's 
theorem. The desired injection energy is then determined 
using an iterative procedure with arrival position as an 
isolation parameter. The methods used to assure the com- 
putation of proper type and class trajectory are described. 

Curves are included showing minimum and constant C 3 
for trajectories from Earth to Jupiter in the 1970-71 syn- 
odic period with flight time, central angle, inclination 
of transfer plane, right ascension of launch asymptote, 
declination of launch asymptote, and arrival energy as 
dependent variables with launch date as independent vari- 
able. A discussion of the modification required to solve 
other interplanetary problems of current interest is also 
included. A detailed comparison of the computer program 
utilizing these procedures with other methods in use indi- 
cates the accuracy is the same, but this method uses 
considerably less computer time. 


SECTION I . INTRODUCTION 

The basic problem presented in this report is the 
time-constrained two-point boundary value problem of calcu- 
lating the trajectory which under the influence of only the 
sun's gravitational force connects the position of the 
launch planet at launch time to the position of the arrival 
planet at arrival time. The methods used to solve this 
problem differ only slightly from those used in earlier 
work [2,3]; however, the way in which the solutions are 
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presented and the way the solutions are used to advantage 
in solving additional problems are unique. 

The planetary orbits are represented by inclined 
elliptical orbits, and a form of Lambert's theorem [9] is 
used to solve the time-constraint. After the conic transfer 
trajectory is determined, the injection energy is calculated 
in terms of the impulsive velocity increment required to get 
from the planetary orbit to the transfer orbit. The object 
of the procedure to be discussed is to obtain a trajectory 
requiring a given injection energy. The position of the 
target planet at arrival is varied as an isolation parameter 
to obtain the desired injection energy. Normally, it is 
desirable to solve the problem over a range of launch dates, 
and it is possible to use the results of each launch date 
as a basis to extrapolate the first guess to the solution 
at the next launch date . 

A brief discussion of the modifications which must be 
made to the procedure in order to solve additional problems 
of current interest is also included. It is felt that the 
main advantages of the procedures presented are speed, high 
versatility, and ability to solve the problem in terms of 
injection energy. Comparisons with some other procedures 
will be made to emphasize these points. 

The description of the problem and the methods used 
to solve it will, as much as possible, be made without the 
use of mathematical equations. However, a complete list 
of the equations used to solve this problem and an indication 
of how they were obtained will be found in the Appendix. 


SECTION II. PLANETARY COORDINATES 

It is desirable to designate the launch date; and 
therefore, it is necessary to have some method of determining 
the planet's positions and velocities as a function of time. 

The most accurate method used involves interpolation between 
values in an ephemeris [3,4]. However, a method which is 
sufficiently accurate for the conic transfer, and a great 
deal faster, involves the use of mean conic elements for 
the planetary orbits [2,4]. 

It is necessary to have five elements for each planet. 

The ones chosen are time at perihelion (t p ), semi -major axis (a), 
eccentricity (e), inclination (i) of launch planet's orbital 
plane to that of the target planet, and the true anomaly of 
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the ascending node (w) of each planet with respect to the 
other. Prom these five basic elements, all additional 
heliocentric constants which may be required can be determined. 


SECTION III. INTERPLANETARY TRAJECTORY 

Once the positions of the launch and target planets 
have been established, the problem is to determine if there 
exists an elliptical trajectory about the sun which will con- 
nect the two points in the required time. There is a one 
parameter family of ellipses which pass through the two points 
with the sun at their focus. Considering only direct or 
counterclockwise flight of less than 360°, at most one of 
these ellipses will have the required flight time. Lambert's 
theorem is used to select the required ellipse. The over-all 
geometry of the planetary transfer is shown in Figure 1. 

Lambert's theorem expresses the transfer time, t, as a 
transcendental function of the radius to the launch and 
target planets (r s and r^), the central angle of flight ($), 
and the semi -major axis of the transfer ellipse (a). A 
complete description of Lambert's theorem and its use is 
given by Breakwell [2]; however, the following description 
should be sufficient to provide a general understanding of 
the theorem. To facilitate the solution of the problem, 
two parameters are defined. The first parameter, K, has 
a range of 0 to 1 and depends only on the planetary geometry 
(i.e., K = f x (rs,ric, $) ). The second parameter, E, has a 
range of -1 to 0 and depends on the semi -major axis of the 
transfer ellipse as well as the planetary geometry (i.e., 

E = f 2 (r s ,rk, $,a) ) . The transfer time, t, is then expressed 
as a function of the two parameters, K and E (i.e., t=f(K, E)). 
The function f is normally expressed in terms of inverse 
trigonometric functions, but for values of E near 0, and t 
small, trigonometric formulation is neither accurate nor the 
most rapid. Thus for E greater than -.2 and t less than t fc 
a series formulation for f is used. The quantity t’ D is 
explained in the following paragraph. 

Plots of t versus E for various values of K are shown 
in Figure 2. The minimum time for an elliptical transfer 
is the time, t a , at E=0. Transit times smaller than t a 
require a hyperbolic transfer. As shown in the figure, for 
values of K greater than 0, there are four values of t for 
each value of E. The following rules are used to determine 
which value is correct for each problem. Branch 1 is used 
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for central angles less than l8o° and Branch 2 for central 
angles greater than l8o°. The time, t^, can be calculated 
by setting E=-l; then the lower portion of the curve Is 
used If t is less than tfc and the upper portion if t is 
greater than t^. Since each portion of the curve is mono- 
tonic, a two-point interpolation scheme can be used by 
providing that interpolated values of E are restricted to 
the range -1 to 0. The first interpolation is made between 
the points (-1, t’ D ) and (E x , t x ) where E x is a first guess 
at the solution and t x is the associated value of t. 
Henceforth, the interpolation is made between the last two 
points calculated until the calculated value of t is within 
the desired tolerance of the actual t. Once this condition 
is reached, the transfer ellipse is completely determined. 


SECTION IV. INJECTION ENERGY 

The injection energy, C 3 , is the energy of the launch 
planet centered conic. It can be calculated as the square 
of the hyperbolic excess velocity vector (v n ). This veloc- 
ity vector is defined, as shown in Figure 3, as the differ- 
ence between the launch planet velocity vector (vp) and the 
velocity on the transfer ellipse at launch time (vt). Two 
implicit assumptions being made here are that the velocity 
addition is impulsive and that it is made at the center of 
the launch planet. The error introduced by these assumptions 
is sufficiently small for the purpose of this problem. 

Graphs of various trajectory parameters versus launch 
date with C 3 as a parameter, such as those in Figures 12 
through 17, have been very useful in mission planning pro- 
jects and in projects involving general surveys of inter- 
planetary trajectories. Previously, these curves were 
generated by running a large number of trajectories at 
discrete intervals of flight time on each of many launch 
dates. Then the constant C 3 can be established using a 
curve fit interpolation scheme. This interpolation must be 
performed on each trajectory parameter on each launch date. 

It seems, however, that it is more efficient and desirable 
to calculate these constant C 3 curves directly. 

It is interesting and necessary to calculate the minimum 
injection energy for a given launch date. The first reason 
is that minimum injection energy could be useful in accom- 
plishing a given mission with maximum payload, or in 
preliminary planning stages, it could give an idea of the 
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absolute minimum energy characteristics required .of the 
launch vehicle. The second reason is that in order to 
calculate constant C 3 values over a range of launch dates , 
it is necessary to know first, for each of the days, that 
the minimum C 3 for that day is less than or equal to the 
constant C 3 value being investigated. A plot of C 3 versus 
arrival planet 's position (©k), for a given launch date, is 
shown in Figure 4. As shown, there are two relative minimums. 
The trajectories for central angles less than l8o° are called 
Type I and those with a central angle greater than l8o° are 
called Type II [3]. Thus, there is a Type I and a Type II 
minimum. 

The method used to calculate the minimum is to take a 
first guess at the 0k which will give the minimum C 3 . Then 
the C 3 corresponding to 0k±A, where a is a small increment 
on 0k, is calculated. Next, a parabola is passed through 
these three points, and the minimum point on this parabola is 
used as the next guess at 0k. A restriction is placed on 0k 
so that it will provide an elliptical trajectory of the pro- 
per type. The two lowest of the last three points plus the 
new point are again fit with a parabola to obtain a better 
guess at 0k, and the procedure is repeated until the minimum 
is obtained to a satisfactory accuracy. Once the solution 
is obtained, it is saved for use in extrapolating the first 
guess on 0k at the next launch date. The minimums for 
successive launch dates are calculated over the required 
range of launch dates . These minimums determine the curves 
shown in Figures 5 through 11. 

As indicated above, the launch date range of constant 
C3 curves is determined from the minimum C 3 curves. For 
example, it might be desirable to calculate a Type I 
curve of 100 km 2 /sec 2 for an Earth -to-Jupiter transit in 
the 1971 synodic period. Referring to Figure 5, Type I 
trajectories with a C 3 =100 km 2 /sec 2 are available only 
between January 11 and February 28. For launch dates out- 
side this range, a higher C3 is required. To illustrate how 
a point on a given C 3 - curve is calculated, it will again be 
necessary to refer to Figure 4 where it can be seen that 
there are two points with a C 3 of 100 km 2 /sec 3 for each type 
trajectory at the indicated launch date. These two solutions 
are” identified as Class I and Class II, Class I always being 
the trajectory with the shorter flight time. It is apparent 
that the two classes are coincident for the minimum C 3 on any 
launch date. The slope of the C 3 curve is used to determine 
which class is which. If the slope is negative, the solution 
is Class I; if the slope is positive, the solution is Class II. 


5 


The method used to isolate a given C 3 value is to take a 
first guess at and to calculate the C 3 for this 8^ and 

for 0k+A* Then a linear interpolation scheme is used 
between the last two points calculated until the desired 
C 3 is established within the given tolerance. During this 
iteration the points are restricted such that the proper 
type and class trajectory will be computed. Again, the 
solution from one case is used to establish the first guess 
at the solution to the following case. Figures 12 through 
17 show typical plots of launch date versus various dependent 
variables for several C 3 values. 

SECTION V. RESULTS AND CONCLUSIONS 

Extensive comparisons between this computer program 
and the JPL and Lockheed programs were made during the 
development and checkout of this deck. The JPL program 
used for comparison was their Heliocentric Conic Program [3] 
which uses ephemerides to determine the positions and veloc- 
ities of the planets. It was found that this deck requires 
about one second to compute each interplanetary trajectory. 

The Lockheed program used for comparison was their Medium- 
Accuracy Interplanetary Transfer Program [2] which uses 
mean conic elements of the planets to determine the planets ' 
positions and velocities. It was found that this program 
requires about .1 second for each interplanetary trajectory. 
Table 1 summarizes the deviation of nine parameters for 100 
arbitrarily chosen trajectories from Earth to Venus, Mars, 
and Jupiter. The only restriction placed on these trajec- 
tories was that the central angle of flight not be between 
170° and 190°. This was done because the rate of change of 
energy and transit plane inclination with respect to plane- 
tary position is very high in the near l8o° central angle 
region. Thus, small discrepancies in planetary position, 
which do appear due to the different methods of calculating 
position, cause large discrepancies in energy and inclination. 
This difficulty is not felt to be a major drawback, first, 
because even though the values are not the same, the trends 
are the same in all three decks. Secondly, the problems occur 
at values of C 3 much higher than those which are practically 
feasible . 

Table 2 shows a comparison between these three conic 
programs and the JPL Space Trajectories Program which is a 
very accurate n-body integration program. The trajectory 
parameters chosen for comparison give an indication of the 
accuracy of the conic programs at both the injection and 
terminal ends of the trajectory. The approximate computation 


6 



times given in the chart indicate the tremendous savings in 
time accomplished by using conic calculations. The time 
per trial given for the integrating deck is for a single 
trajectory calculated with a given set of initial conditions. 
Using initial conditions from one of the conic decks, the 
integrating deck will isolate the desired terminal condition 
in about 20 trials. If less accurate initial conditions are 
used, more trials will be required. 

The graphs of Figures 5 through 17 are examples of the 
type of material that can be produced directly from the out- 
put of this program. Each of the graphs with constant C 3 
curves was drawn using 450 points with a total computation 
time on the IBM 7094 of about 31 seconds. The minimum C 3 
curves used an additional 84 points and required an additional 
16 seconds of computation time. It is estimated that it would 
require at least 1000 trajectories to produce these same curves 
by the old method of interpolating between points in a matrix 
of values. Additional computations would also be required to 
perform the interpolation for constant values of c 3 . It is 
apparent then that not only is the program more efficient at 
solving the basic interplanetary transfer problem, but, 
additionally, it is more efficient at solving the over-all 
problem of determining constant C 3 curves because of the 
new approach to the problem. 

SECTION VI. OTHER USES 

With minimum modification, the principles used in this 
procedure can be used to solve other interplanetary trajectory 
problems of interest. For interplanetary orbiting or landing 
missions, total velocity increment, that is, injection veloc- 
ity plus braking velocity at the target planet, is of more 
interest than C 3 . Thus, by calculating this total velocity 
and using it in place of C 3 , the same procedure can be used 
to generate minimum and constant total velocity data. In 
much the same way, payload mass at the target planet can be 
approximated by impulsive calculations if the characteristics 
of the vehicle system are known. 

Round -trip fly-by missions and "grand tour" missions, 
where no major propulsion is used after injection, can also 
be handled quite easily. For example, suppose it is desirable 
to fly a round-trip fly-by mission to Mars. The first step is 
to calculate minimum energy requirements for Earth-t o -Mars 
and for Mars -to-Earth trajectories. Then the constant C 3 range 
for launch at Earth can be chosen and a trajectory from Earth- 
to-Mars with the desired C 3 can be established. Now the energy 
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of the Mars centered conic at arrival can be determined, 
and this can be used to calculate a Mars -to -Earth leg with 
the same energy, provided this energy is not less . than, the 
minimum energy for this particular Mars passage date. The 
hyperbolic excess velocity vector of approach to Mars is 
then rotated into the hyperbolic excess velocity vector of 
departure by calculating the proper radius vector from the 
center of Mars to the pericenter of the Mars centered conic. 
This same procedure could be used for "grand tour" type 
trajectories except that more than two legs would be involved, 
and the conics would have to be matched at each of the planet 
passages . 

It is felt that the speed and flexibility of this program 
will make it very useful in approximating trajectories for 
any interplanetary mission, or, in a broader sense, in cal- 
culating free flight transfer trajectories between any two 
orbits around any central body. 
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COMPARISON BETWEEN N-BODY AND CONIC CALCULATIONS 
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Nodal Line 


FIG. 1. PLANETARY TRANSFER GEOMETRY 
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(km 2 /$ec 2 ) 



INJECTION ENERGY VERSUS ARRIVAL PLANET'S POSITION 
TO JUPITER TRAJECTORIES FOR LAUNCH ON JAN. 31, 1971 










(joas/jun) 



17 


TG. 5. LAUNCH ENERGY VERSUS LAUNCH DATE 
MINIMUM ENERGY EARTH TO JUPITER TRANSFER 








FIG. 6. FLIGHT TIME VERSUS LAUNCH DATE 
MINIMUM ENERGY EARTH TO JUPITER TRANSFER 
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16. 7. CENTRAL ANGLE VERSUS LAUNCH DATE 
MINIMUM ENERGY EARTH TO JUPITER TRANSFER 









FIG. 8. INCLINATION OF TRANSFER ORBIT VERSUS LAUNCH DATE 
MINIMUM ENERGY EARTH TO JUPITER TRANSFER 




RLA (deg) 
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FIG. 9. RIGHT ASCENSION OF LAUNCH ASYMPTOTE VERSUS LAUNCH DATE 
MINIMUM ENERGY EARTH TO JUPITER TRANSFER 








Launch Date 



DECLINATION OF LAUNCH ASYMPTOTE VERSUS LAUNCH DATE 
MINIMUM ENERGY EARTH TO JUPITER TRANSFER 





(10 4 km 2 /sec 2 ) 
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FIG. 11. ARRIVAL ENERGY VERSUS LAUNCH DATE 
MINIMUM ENERGY EARTH TO JUPITER TRANSFER 
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CENTRAL ANGLE VERSUS LAUNCH DATE 
iARTH TO JUPITER TRANSFER 























RLA (deg) 
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ASCENSION OF LAUNCH ASYMPTOTE VERSUS LAUNCH DATE 
EARTH TO JUPITER TRANSFER 
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DECLINATION OF LAUNCH ASYMPTOTE VERSUS LAUNCH DATE 
EARTH TO JUPITER TRANSFER 
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ARRIVAL ENERGY VERSUS LAUNCH DATE 
»RTH TO JUPITER TRANSFER 













FOR CONSTANT AND MINIMUM C 
































APPENDIX 


This appendix is designed to give a mathematical treat- 
ment to the material presented in the main report. It is not 
designed as a complete description of the computer program. 

A userls guide to be published by the MSFC, Computation 
Laboratory will satisfy this requirement. The simplified 
flow diagram of Figure A1 shows the order used to perform 
the required calculations. The calculations performed in 
each block are listed in the descriptions of the subroutines 
which follow. 


PLANET POSITION I 

Using date and the mean conic elements of the planets , 
this routine calculates the position of the desired planet 
and the angular position at which the other planet will be 
at the transition between Type I and Type II trajectories [6]. 

M j = Hj^j-tpj) 0 < M < 2tt 

Ej is calculated by iteration on the equation 

0 < M < 2tt 

0 _< 0 . < 2-rr 
d 

0 a- < 2rr 
0 < 9 X < 2tt 


M, = 


e j 


Ej-ej sin Ej 


1+e 


2 arc tan tan Ej/2 

J 


r j 


“j 


A j ( 1 -e j cos E j ) 

e J - 


tan 


e„ = w, . + arctan l •— 0 . 
x k,j Vcos i 






PLANET POSITION II 

Using the true anomaly of the desired planet, the mean 
conic elements of the planets, and the date and position of 
the other planet, this routine calculates the radius and date 
of the desired planet, the angular separation of the two 
planets, and the flight time [6]. 
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DT = T k - Tj 

LAMBERT CALCULATION 

Using the positions of the two planets and flight time., 
this routine calculates the required transfer ellipse [2,6]. 

The chord from the launch to arrival point is 



s = ( r j + r k + °)/ 2 


E — -1 < E < 0 

2a' ~ “ 

where a' is a first guess at the semi -major axis of the 
transfer ellipse. 

C 

K = 1 

s 

Normalized time for a parabolic transfer is 
t a = 1(1 - P K*) 

where P = +1 if $ < 180° and P - -1 if $ > l80°. 
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Normalized required transfer time is 


86400 x > f2jl 
t -s- ra 

and time for a minimum heliocentric energy transfer is 

TT 

t b P(arcsin W - %(l-K) ). 

2 

If t < t b , Q = +1; if t > t b , Q = -1. Then the desired 
E is found by iterating the solution on one of the following 
transcendental functions 

t' = (-E)"»(|* + Q(“ + arcsin - ^E(l+E) ) 

- P(arcsin 1-E K - A|e K(l+E K) ) j . 

If Q > 0 and E > -.2, a series solution, obtained by 
expanding each term in a series and adding term by term, 
is used. 

N 

t' = c i ( x “ p K^ +1 )E 1_1 . 
i=l 

The first ten coefficients are 

C x = . 66666 667 

C 3 = .10714286 
C 5 = .049715909 
c 7 = .030078125 
C 9 = .020671644 

The next guess at E is found by making a linear interpolation 
between (-l,t b ) and (E,t'). Then a new value of t ! is cal- 
culated and additional linear interpolations are made between 
the last two calculated points until t' is within the desired 
tolerance or t. 


C = - 2 

C l = -!o69444444 

c 6 = -.037860577 

C 8 = -.024643842 

Cio= -.017663865. 
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Then the elements of the transfer ellipse are calculated 
as follows 



arctan 


(l+E)(s-r k >/(l+EK)(s-r j ) 



TT 

2 






. -1 

i x = sm 

x = V 1 cos Y 

y = V x sin 1 1 sin y 
z = V 1 cos l x sin y 


sln h.k sln g k 

sin $ 




0 < 9 < 2rr 


- < 1 < - 

2 “ “ 2 
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C 3 ROUTINE 


Using the velocity of the probe and the mean conic 
elements of the planets, this routine calculates the velocity 
of the planet, the velocity of the probe with respect to the 
planet, and the planetocentric energy of the probe [6]. 

_ <TPP,i 



sin e 
PP— 


* , 

x' = x - r. 

d 

y' = y - 0j 

z 1 = z 

C 3 = x ' 2 + y' 2 - + z 12 


INTERPOLATION FOR CONSTANT C 3 

This is not a subroutine but is incorporated in the main 
line of the program. The slope of C 3 vs is checked to be 

sure the interpolation is being made for the proper class 
trajectory. Then a linear interpolation between the last 
two points is made to determine the which will produce the 

desired C 3 . If both values of C 3 used for the interpolation 

are above the desired value., the calculated change in 9^ is 

increased slightly to avoid creeping on the solution. 

INTERPOLATION FOR MINIMUM C 3 

This routine uses three points on the C 3 vs e curve to 

interpolate for the 0- which will produce the minimum point 
on the curve [7]. 

The three points used for interpolation are designated 
fey (e k , C 3x ), (e^, C 3a ) and (e ka , c 3a ), and are arranged 
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such that C 3 <C 3 < Co . Then to reduce scaling problems 

a transformation of coordinates is made such that the origin 


is at the first point. 



X 

II 

CD 

X 

M 

1 

CD 

X 
(— 1 

II 

o 

yi = °3 X 

O 

II 

<5* 

o 

X 

to 

II 

CD 

1 

CD 

FT 

y 3 = c 3s 

- c 3 

3 x 

X 

CO 

II 

CD 

x- 

CO 

1 

CD 

x 

H* 

y 3 - c 33 

1 

o 

CO 

H 


The equation of a parabola through these three points is 


x(x-x 3 ) x(x-x 2 ) 

y = — : r Vs + 


x a (x 2 -x 3 ) 


x 3 (x 3 -x a ) 


y 3 




2x-x c 


2x-x ; 


dx -x 3 (x a -x 3 ) + x 3 (x 3 -x 2 ") ya °' 

Thus the minimum is at 


x = 


x 8 2 y 3 - x 3 s y 8 
2(x 8 y 3 - x 3 y 2 } 


and 


EXTRAPOLATION ROUTINE USING LAGRANGE FORMULAS 

Using the last M solutions to the problem, this routine 
extrapolates the first guess on e k and a for the solution 

to the problem at the next launch date [7]. 

The values of Tj , a, e^. are saved in a table and 
designated as T 1 , a i , . 

For the second point Tj is incremented and a is set 
equal to a x and is set equal to . 
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If two or more points are in the table is incremented 

and the following formulas are used to calculate first guesses 
for a and 0 k , 

M 


N = (T - T.) 


i=l 
M 


D „ = 


= TT - T ±) 1 ¥• t i = l, M 


1=1 


c„ = 


N 


* W? 


M 


'k 


- L °e V 


1=1 

M 


= Y~ C a 
Z & i 


a=l 


A = 1, M 


VELOCITY CALCULATION ROUTINE 


Using the elliptical elements of the transfer trajectory, 
the mean conic elements of the planets and the position of 
the arrival planet, this routine calculates (l) velocity of 
the probe at arrival, (2) velocity of arrival planet at 
arrival, (3) inclination of transfer conic to arrival planet 
conic, (4) velocity of probe with respect to arrival planet 
and (5) energy of the planetocentric conic [6]. 


e k = 


PPk 

r k 


e k sin 9 k 
r k = PP£ 

i 2 - sin i, . sin a . /sin $ 

k,j j- 

0 2 = PP/ r k 

r 2 = VpT e(sin ej/pp 
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- r k 


V x = 

Vy = 02 COS 1 3 — 9^. 
V 3 = 0 3 sin i 2 
C 3 A = V x 2 + Vy 2 + V z 2 


SUBROUTINE ROTATE 

This routine rotates a vector in a coordinate system 
where the x-axis is along r, the z-axis is perpendicular to 
the planet's orbital plane, and the y-axis completes a right- 
handed system into a coordinate system analogous to the 
earth's "Inertial Cartesian Equatorial System." That is, a 
coordinate system where the x-axis is in the direction of 
the "vernal equinox" (sun ascending node of the planet's 
equatorial plane), the z-axis is perpendicular to the planet's 
equatorial plane, and the y-axis completes a right handed 
system. After the rotation, the routine calculates the right 
ascension and declination for an "Inertial Spherical Equatorial 
System" [2], 


X * 


/cos 

( 0 i- T i 

) 

-sin(e 1 -T i 

) 

0 ^ 

y T 

= 

sin 

( e i- T i' 

)cos eij_ 

COs(e i -T i 

)cos eij_ 

-sin eij_ 

z 1 


i^sin 

b- 

1 

•H 

CD 

)sin eij_ 

COS(e i -T 1 

)cos ei^ 

cos el^J 


RLA = 

arctan 

(y'/x' ) 





(raa) 







DLA = 

arctan 

. ( Z ' /Vx ' 2 + 

y 12 ) 
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